from __future__ import division
import numpy as np
from Constant_Physics import Physics
import src_fortran.froutines as fr

###
# Check matrix
def Check_population_conservation(input_matrix,matrix_name,print_population):
    """
    Input:
        - matrix has to be square matrix
    Output:
        -Binary: True if conservative, else False
    """
    # First, try the self.NonResonantRates matrix
    # For each column, the sum of all the elements should be zero
    conservative_matrix=True
    dim=int(np.sqrt(input_matrix.size))-1
    nrr=np.zeros([dim])
    for k in range(dim):
        nrr[k]=np.sum(input_matrix[:,k])
        if np.abs(nrr[k]) >= 1e-8:
            conservative_matrix=False
    if conservative_matrix==False:
        #print "\nThe ",matrix_name," matrix is not conservative:"
        if print_population==True:
            print nrr
        return False
    else:
        #print "\nThe ",matrix_name," matrix is conservative."
        if print_population==True:
            print nrr
        return True


###
# Obsolete
class Lorentzian():
    """
    Computes Lorentzian
    """

    def __init__(self,derivative=0):
        self.derivative=derivative

    def Eval(self, coeffs, evalpts):
        """
        evaluate lorentzian with given coeffs over given evalpts 
        coeffs = (A0,E0,G0,n)
        """
        A0, E0, G0 = coeffs
        x = np.asarray(evalpts) #XXX: requires a numpy.array 
        if self.derivative==0:
            return A0 * (G0 / (2 * np.pi)) / ((x - E0) * (x - E0) + (G0 / 2) * (G0 / 2))
        elif self.derivative==1:
            return (16*A0*G0*(E0-x))/(np.pi*(G0**2+4*(E0-x)*(E0-x))*(G0**2+4*(E0-x)*(E0-x)))
        elif self.derivative==2:
            return -(16*A0*G0*(G0**2-12*(E0-x)*(E0-x)))/(np.pi*(G0**2+4*(E0-x)*(E0-x))*(G0**2+4*(E0-x)*(E0-x))*(G0**2+4*(E0-x)*(E0-x)))



if __name__ == "__main__":
    lorentzian=Lorentzian(0)
    intensity=1
    center=9076
    width=3
    coefs=intensity,center,width
    print lorentzian.Eval(coefs,9076)
